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Abstract 


Full  implicit  finite  differencing  was  used  to 
approximate  transient  heat  conduction  in  two  spacial 
dimensions.  The  resulting  set  of  linear  simultaneous 
equations  was  solved  using  successive  over-relaxation 
iterative  methods.  Errors  between  the  exact  matric 
solution  and  the  iterated  solution  were  computed  and 
compared  with  several  methods  of  determining  error 
approximations  that  use  the  displacement  between  iterations 
and  an  associated  parameter.  Fundamental  parameters  of 
the  iteration  matrix,  such  as  the  spectral  radius,  'were 
tested  as  parameters  in  the  error  approximation  methods. 

The  error  methods  were  compared  'with  respect  to  accuracy, 
ease  of  computing,  and  computer  resources  required. 

The  nodal  density  used  in  the  numerical  approximation 
to  the  physical  problem  was  compared  with  the  spectral 
radius  of  the  resulting  iterative  matrix  problem.  The 
spectral  radius  increases  as  nodal  densit;'  increases. 


I  Introduction 


Background 

In  practice,  few  physical  problems  involving  partial 
differential  equations  can  be  solved  analytically.  Most 
initial  and  boundary  conditions  encountered  in  practice 
lack  the  required  symmetry.  Alternatively,  numerical 
approximations  can  be  used  to  generate  rnatric  equations 
that  have  solutions  approximating  the  physical  solution. 

By  increasing  the  nodal  density  used  in  the  numerical 
approximation,  the  error  between  the  rnatric  and  physical 
solutions  is  decreased,  but  at  the  expense  of  increasing 
the  number  of  rnatric  elements. 

Iterative  techniques  are  normally  used  to  solve  rnatric 
equations  because  iterative  techniques  are  less  susceptable 
to  error  instabilities,  they  will  not  accumulate  truncation 
errors,  and  they  require  fewer  computer  operations, 
especially  if  many  rnatric  elements  equal  zero. 

Many  mathematicians  are  troubled  by  the  fact  that 
they  do  not  know  when  to  stop  the  iterative  process.  The 
error  between  iterated  and  rnatric  solutions  is  normally 
unknown. 

Problem 

The  primary  purpose  of  this  thesis  is  to  compare 
methods  of  determining  error  limits  and  approximations  with 
respect  to  accuracy,  computational  ease,  and  required 
computer  resources. 

The  secondary  purpose  is  to  determine  what  effect 
increasing  the  nodal  density  has  on  the  iterative  conver¬ 
gence  rate. 

Scope 

This  thesis  is  limited  to  the  study  of  transient  heat 
conduction  in  two  spacial  dimensions.  Only  the  full 
implicit  numerical  approximation  to  the  physical  problem 
was  used.  The  primary  iterative  method  examined  was  the 


successive  over-relaxation  method 


j 

Approach  and  Presentation 

Section  II  discusses  in  detail  the  theory  that  is 
necessary  to  attempt  a  solution  to  the  problem  of  iterative 
error  analysis.  First,  the  matric  approximation  to  the 
physical  problem  is  explained,  followed  by  the  solution  of 
implicit  matric  equations  by  iterative  techniques.  The 
spectral  radius  is  covered  in  preparation  for  examining  the 
convergence  properties  of  iterative  matric  equations.  The 
concept  of  norms  is  introduced  to  help  quantitatively 
define  iterative  error  limits. 

Section  III  derives  the  sample  problem  investigated 
and  then  the  computer  codes  are  examined  which  produce  error 
related  data. 

Section  IV  presents  the  numerical  data  in  graphical 
form  and  trends  are  discussed.  Both  spectral  radius 
versus  nodal  density  and  iterative  error  approximations  are 
discussed. 

Section  V  draws  conclusions  and  gives  recommendations 
based  on  the  data. 
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II  Theory 

■*,  W 

Matric  Approximation  to  Physical  Problem  (5:55-56,59, 

106-109 ; 9 ) 

In  this  section  we  discuss  the  construction  of  matric 
equations  that  approximate  transient  heat  conduction.  The 
two  dimensional  transient  heat  conduction  problem  can  be 
formulated  into  the  following  parabolic  partial  differ¬ 
ential  equation. 

/tA1  *  (l/K)  )t/ 

In  many  instances,  this  equation  with  its  associated 
boundary  and  initial  conditions  cannot  be  solved  analyti¬ 
cally.  One  method  of  simplifying  this  continuous  problem, 
is  to  make  the  equation  and  the  domain  of  the  variables 
discrete.  First,  a  grid  is  superimposed  over  the  variable 
domain.  The  function,  T,  is  then  expanded  in  a  Taylor 
series  about  selected  adjacent  grid  points.  Linear  combi¬ 
nations  *of  these  are  used  to  form  difference  approximations 
to  the  differentials  (5:59).  The  following  are  six  standard 
difference  operators,  the  differentials  they  approximate, 
and  the  associated  truncation  error  (5:55-56). 

First  forward  difference 

AT(x)  -  T  (*  ♦/))  -  T(h) 

)t/)x  =  AT  A  ♦ 

First  central  difference 

$  T60  =  T (x*h)  -  T(x-h ) 

>  zr/ix  =  $>  t/ 2h  i-  <y(^) 
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First  backward  difference 


7tm  = 

iT/ix  =  VT/h  *  o(h ) 

Second  forward  difference 

A*  17*)  =  -  2T(x+h)  t  j  60 

=  a"  rA7  ♦ 

Second  central  difference 

S*  =  r  -ztW  <•  rq-O 

ilT/)x‘  =  fT/i,'  *  oft1) 

Second  backward  difference 

7*TW  =  T60  -ZT(*-h)  +  T(*-Zh) 

y  t/jx1  =  v'T/y  +  <>(7,) 

where  h  is  the  grid  interval  in  the  x  coordinate. 

Replacing  differentials  with  finite  differences,  a 
set  of  linear  difference  equations  is  created.  One 
equation  exists  for  each  spacial  grid  point,  or  node, 
interior  to  the  boundary.  These  equations  are  then  written 
compactly  as  a  matric  equation. 

U  -  A  V  +  b 


whore 


A  =  matrix  of  coefficients 
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W  =  vector  of  knovm  (or  unlcnovm)  nodal  values  of 
the  function  T 

V  =  vector  of  unknown  (or  known)  nodal  values  of 
the  function  T 

Id  =  vector  of  known  constant  elements  derived 
from  boundary  conditions. 

This  matric  equation  is  solved  once  for  each  discrete  time 
step. 

Selection  of  a  difference  operator  to  approximate  a 
given  differential  is  governed  by  truncation  error,  error 
stability  and  computational  ease  of  the  resulting  matrix 
equation.  Several  combinations  of  difference  approximations 
exist  that  transform  the  parabolic  partial  differential 
equation  into  two  basic  forms  of  matric  equation.  One  is 
the  explicit  form  and  the  other  is  the  full  implicit  form. 

Each  difference  equation  of  the  explicit  form  expresses 
a  single  unknown  nodal  quantity,  T,  in  terms  of  several 
knovm  nodal  quantities.  In  the  matric  equation 

u  =  A  v  +  b 

u  becomes  the  unknown  vector  and  v  becomes  the  known  vector. 
Solution  requires  simple  matric  multiplication.  Error 
stability  is  conditional,  depending  on  the  time  and  space 
grid  intervals. 

Each  difference  equation  of  the  full  implicit  form 
expresses  several  unknown  nodal  quantities  in  terms  of  a 
single  knovm  nodal  quantity.  In  the  matric  equation 

u  -  A  v  +  b 

* 

*  i, 

u  now  becomes  the  knovm  vector  and  v  becomes  the  unknown 
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vector.  Solution  requires  a  method  equivalent  to  matric 
inversion.  Error  stability  is  unconditional. 

A  weighted  combination  of  the  implicit  and  explicit 
forms  (5:108-109)  can  be  used  to  reduce  truncation  error. 
The  special  difference  operators  of  the  implicit  form  are 
weighted  by  a  factor  <x  ,  v/hile  the  spacial  difference 
operators  of  the  explicit  form  are  weighted  by  a  factor 
I  -  x  .  The  combination  form  is  conditionally  stable  if 

OM<  1/2 

and  unconditionally  stable  if 

1/2  *  *  £  I 

When  is  set  equal  to  H,  the  combination  form  is  called 
the  Crank-IIicholson  method.  Crandall  found  the  *><  that 
minimizes  the  truncation  error.  Truncation  error  can 
also  be  reduced  by  decreasing  the  spacial  grid  interval. 

Solution  of  the  Implicit  matrix  Equation 

Several  methods  exist  for  solving  the  implicit  matric 
equation 


A  x  =  b 

directly.  The  matrix  can  be  inverted  (4:125-129)  to  form 
the  explicit  matrix  equation 

x  =  A'’  b 

the  vector,  x,  can  be  found  by  using  Gauss  reduction 
(5:9-11),  or  if  the  matrix  is  block  tridiagonal,  the  Thomas 
method  of  matrix  factorization  (5:46-49)  can  be  used  to  find 
x.  The  Thomas  method  is  a  special  case  of  Gauss  reduction. 
Unfortunately,  these  direct  methods  suffer  from  a 
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serious  flaw.  The  number  of  algebraic  operations  required 
is  proportional  to  a  power  of  the  number  of  elements  in  the 
unknown  vector,  x.  Each  algebraic  operation  produces 
truncation  error  which  accumulates  as  the  number  of 
operations  increase.  When  the  number  of  unknown  vector 
elements  is  large,  the  accumulated  error  can  be  significant 
(5:111) . 

Alternately,  the  implicit  matric  equation 

A  x  =  b 

can  be  rewritten  as  (12:236) 

x  -  G  x  +  c 

An  iterated  solution  can  then  be  attempted 

(n+i)  -  /-  (*)  ^  s- 

X  "  tJ  x  +  C 


where 


(r>) 

X  =  nth  iterated  vector 

O  =  iteration  matrix 
C  =  known  constant  vector 

The  advantage  here  is  that  truncation  error  accummulated 
during'one  iteration,  is  not  passed  on  to  the  next.  Also, 
by  using  an  equivalent  algorithm,  the  number  of  algebraic 
operations  are  reduced  if  some  elements  of  matrix  A  equal 
zero  (5:111) . 

V/e  now  define  the  displacement  vector  as  the  dif¬ 
ference  between  two  successive  iterated  vectors. 


s' 
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J<">  =  -  X(n) 

The  displacement  will  be  used  here  to  describe  the  iterative 
algorithms . 

The  three  iterative  methods  most  commonly  used  are  the 
Jacobi,  Gauss-Seidel,  and  successive  over-relaxation  (SOR) 
methods.  Their  algorithms  are  derived  by  rewritting  the 
implicit  matrix  in  terms  of  its  strictly  upper  triangular, 
diagonal,  and  strictly  lower  triangular  matrices.  They  are 
denoted  by  U,  D,  and  L  respectively  (12:234). 

Ax  -  b 

A  =  U  +  D+L 
(U  +  D*L)  x  =  b 
L^x+Dx  +  Lx  =  b 
Dx  =  -  U*  -  L  x  +  b 

The  Jacobi  method  expresses  the  (n+l)th  iterative  vector 
elements  exclusively  in  terms  of  the  nth  iterative  vector 
elements  (12:228). 

Dx'"-l)  * -Uxcn'  -  L  x("‘  *  b 
fix'"'0  -  Dx<n)  *  -  Vxw-Dxw-Lx'’'  *b 
x""0  -  X<0>  =D"(-  Ux^-Dx’-’-Lx^'b) 
J(n>  =  D"  (- U  Lx^’  +  b) 

The  Gauss-Seidel  method  attempts  to  accelerate  convergence 
by  using  the  (n+l)th  iterative  vector  elements  as  soon  os 
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they  become  available  (10:229) 


D xM  =  -UxW-j.xM*  b 
D  X  -  D  xCn)  =  -  Ux‘"’- DxM-  L  x1’'"*  +  b 
X^'°  -  X  M  --  D“  (-U  X  -  D  X  -  L  X  *  b  ) 

dfn)=  D"  C-  D  x  fo) -  L  x f*'°  +  b) 

The  SOR  method  attempts  to  further  accelerate  convergence 
by  multiplying  the  Gauss-Seidel  displacement  by  an  accel¬ 
eration  constant,  v/  (10:230). 

J'"'-  w  D"  (-  Ux"')- D  x<n>-  L  <•  b) 

The  iterative  matric  equations 

fn  +  ij  ^/j) 

X  =  C?  x  *  c 

that  are  equivalent  to  the  above  algorithms  are 

Dx'"'0-  -Dx'”-  L  b 

D  x  =  (-U-Ox°*  *  b 
x^'>=  D"(-U-l)km  *  D"b 

for  the  Jacobi  method, 

D  x1"0  =  -Uxm-L  x(n,°  *  b 
D*w*  L  x^*°  -  -Ux^'b 
(D  +  Ox(~°  =  -lixw  »  b 
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x*"'  ‘  (D^O"  (-U)  x'"  >  (d<l)"  fa 

for  the  Gauss-Seidel  method,  and 

XH-XW=  **  D"(-Ux<n,-D*M-L 
Dx""1-Dx"”=-w  Uxw-^Dxw-»L  x^*°^b 

*Dxw-  wD*wt»b 

(D->»-Ox<'’'',=  (-~U  +  D-  ~D)  xf”>  +  >vfc 
x  ^  “  (  D  r  **  L  )  w  U  r  D~  w/t)  x  ^  h  ^D^tvO  wb 

for  the  SOR  method. 

Spectral  Radius 

One  must  understand  the  concept  of  matric  spectral 
radius  in  order  to  study  the  properties  of  iterative  matric 
solutions.  The  spectral  radius  will  be  defined  first. 

Then  a  generalized  method  for  finding  the  approximate 
spectral  radius  and  a  method  for  determining  the  accuracy 
of  that  approximation  will  be  discussed. 

The  concept  of  spectral  radius  is  closely  related  to 
the  concepts  of  eigenvalue  and  eigenvector.  If  a  matrix 
operates  on  a  vector  and  results  in  a  scalar  constant  times 
the  original  vector,  then  that  scalar  constant  is  an  eigen¬ 
value  of  the  matrix  and  the  vector  is  an  eigenvector 
associated  with  that  eigenvalue.  The  equation  is  an  eigen- 
equation. 


M  x  *  A  x 

Unless  a  matrix  is  square,  matrix  multiplication  will 
cause  the  resultant  vector  to  have  a  dimension  different 
from  the  original  vector.  Therefore  eigenvalue  and  eigen- 
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vectors  are  associated  with  square  matrices  only.  The 
number  of  eigenvalues  associated  with  a  square  matrix  is 
equal  to  the  dimension  of  the  matrix.  A  three  by  three 
matrix  will  have  three  eigenvalues. 

If  an  eigenvalue  differs  from  all  other  eigenvalues, 
it  is  termed  distinct.  If  an  eigenvalue  is  equal  to  some 
other  eigenvalue,  it  is  termed  multiple.  The  number  of 
eigenvalues  with  equal  value  is  termed  the  eigenvalue 
multiplicity. 

All  eigenvectors  associated  with  a  given  eigenvalue 
form  a  vector  subspace.  The  number  of  linearly  independent 
vectors  needed  to  span  that  subspace  is  no  larger  than  the 
multiplicity  of  the  associated  eigenvalue  ( 5 : 23-31 ; 1 : 168-171 ) 
Associated  eigenvectors  may  not  be  able  to  span  the  vector 
space  for  non-symmetric  or  non-hermitian  matrices.  Ad¬ 
ditional  principal  vectors  may  have  to  be  added.  A 
distinct  eigenvalue  has  a  multiplicity  equal  to  one.  All 
eigenvectors  associated  with  it  can  be  spanned  by  one  eigen¬ 
vector.  The  eigenvectors  are  all  scalar  multiples  of  each 
other.  (3:131) 

Mx=>x 

a  M  x  -  a  A  x 

Max  -  X  a  x 

My  =  X  y 

Two  of  these  eigenvectors  have  unit  length.  One  is  the 
negative  of  the  other.  (9:35) 

An  eigenvalue  having  a  multiplicity  equal  to  three 
might  need  one,  two,  or  three  linearly  independent  vectors 
to  span  the  eigenvector  subspace  associated  with  it.  Any 
linear  combination  of  eigenvectors  associated  with  the  same 
eigenvalue  is  itself  an  eigenvector  (1:169). 
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ii#  ✓*■»**-« 


ai  M  x  ♦  b  My  =  a  X  x  t  b  X  y 


Max  +  M  \o  y  -  X  a  x  t  A  by 

M  (ax  +  by)  -  X  (dx  *■  by) 

If  an  eigenvalue  is  not  real,  then  its  complex  con¬ 
jugate  is  one  of  the  other  eigenvalues.  The  associated 
eigenvector  subspaces  are  complex  conjugates  of  each  other 
(5: 23-24) . 


M  x  =  X  x. 

M  x*  -  X*  x* 

Zero  is  a  perfectly  acceptable  eigenvalue  that  is 
associated  with  not  all  matrices.  It  will  have  an  associated 
eigenvector  subspace  distinct  from  nonzero  eigenvalues. 

The  matric  spectral  radius  can  now  be  defined  as  the 
modulus  of  the  maximum  modulus  eigenvalue, 

0  =  1x4  *lxd  *  *1X4 

where 

p  =  spectral  radius 


•  =  matric  eigenvalues 

=  dimension  of  square  matrix 


Several  methods  exist  that  can  be  used  to  find  the 
eigenvalues  and  thus  the  spectral  radius  of  a  given  matrix. 
Each  method  has  advantages  and  disadvantages.  One  method  is 
to  solve  the  characteristic  equation  associated  with  the 
matrix. 


|  M-  I  X|  -  O 

V/hen  formulating  this  equation  into  standard  polynomial 

form,  truncation  error  associated  with  calculating  the 

coefficients  becomes  significant  as  order  is  increased 

(3:130).  n  factorial  terms  must  be  summed  to  evaluate  a 

determinant  of  order  n  (11:622).  For  example,  a  deter- 
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minant  of  order  20  would  require  the  summation  of  2x10 
terms!  If  the  resulting  equation  is  quartic,  cubic, 
quadratic,  or  linear,  it  may  be  solved  analytically  (3:104- 
107).  If  it  is  of  higher  order,  one  must  resort  to  itera¬ 
tive  techniques. 

The  eigenvalues  of  some  matrices  can  be  found  by 
inspection.  The  eigenvalues  of  diagonal  matrices,  are  the 
diagonal  elements,  ( 5 : 26 ; 1 : 171-172 ;3 : 13 2) 

The  power  method  is  used  to  generate  a  sequence  of 
iterated  vectors  (5:140)  from  which  several  schemes  have 
been  devised  to  extract  the  spectral  radius  and  an  asso¬ 
ciated  eigenvector.  ( 7 : 5 . 1-5 . 8 ; 5: 144 ; 10: 460)  The  power 
method  consists  of  multiplying  an  initial  trial  vector  by 
the  matrix.  The  result  is  multiplied  by  the  matrix  again 
and  so  on.  Essentially  the  initial  vector  is  multiplied 
by  an  integer  power  of  the  matrix. 

MX.  =  X, 

M  x,  -  x  t 
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M  x„„  -  x„ 

M"  x.  -  xM 

The  following  discussion  will  be  limited  to  real 
symmetric  matrices  although  the  results  can  be  generalized 
to  Hermitian  and  other  special  rnatric  catagories.  All 
eigenvalues  associated  with  real  symmetric  matrices  are 
real.  A  set  of  n  orthogonal  eigenvectors  exists  that  spans 
the  vector  space  associated  v/ith  the  n  by  n  matrix.  (5:24- 
25 ; 1 : 168-169) 

The  spectral  radius  of  a  matrix  is  estimated  by  using 
the  generalized  quotient  method.  An  initial  vector  is 
expanded  in  terms  of  the  eigenvectors. 

x  -  c,  e,  ^  +  + 


where 


C-  =  component  coefficients 

G;  =  eigenvector 

n  =  rnatric  dimension 

^  =  initial  trial  vector 

The  kth  iterated  vector  is  formed. 

x*’  x*” 

*  Mk  i  (c;c;) 

=  |  (c:MkC:) 

I  -  I  ' 
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Me  =  X 


(c:  \he.) 


The  eigenvalues  can  be  ordered  such  that 


lx, I  *  |  Xj  * 


l\J 


Assume  that  the  maximum  modulus  eigenvalue  is  of  multi¬ 
plicity  5^1  .  Only  if  the  negative  of  the  maximum 

modulus  eigenvalue  is  not  an  eigenvalue  itself 


-  A,  *  A 


for  all  i 


can  the  spectral  radius  be  approximated  by  forming  a 
quotient  using  two  vectors  that  are  separated  by  one 
iteration. 

*ft*°*A  r[i(c;e;)'£  (C;[  X:/xf"c;)] 

I -I  l  =  S«l  1 

x(k)=A,k[I  (cie;)>Z  (c,  [A;/ A,]ke;)J 

i  - 1  t  =  5* i 


Since 


then 


|  X;  /  A,  I  *  I  for  I  >  5 


fr-[X,./ A,]1”'  -  kT "  [  x.  /a,]1"  -  o 


and  as  k  increases  toward  infinity 


A,  =  vTx<l'")/(  vT  x(k))  -  ^ 


where  v  is  an  arbitrary  vector.  In  the  limit 
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\  /•'Ti  f  t  (k*i)  / / 

X,  *  k—  L  v  x  /( 


t  (k) 
V  X 


)] 


Assume  that  an  eigenvector  of  multiplicity  f-O 
is  equal  to  the  negative  of  the  maximum  modulus  eigen¬ 
vector. 


-X,  =  X; 


for 


i  =  5  ♦  I  /  S  +  r 


The  spectral  radius  can  always  be  approximated  by  forming 
a  quotient  usinq  two  vectors  that  are  separated  by  two 
iterations 


\r  [  1  (c-.e.t  +  (- ')“‘f  (c:e;) 

i  z  f  i  -5  +■  i 


because 


(.,)*'*  =  (->)«  *  (- 1 /" 

(7: 5. 1-5. 8 ; 10: 460-462; 5: 144) 

In  normal  practice,  the  kth  iterated  vector  is  scaled 
in  some  way  before  iterating  to  find  the  next  iterated 
vector.  In  this  way  working  with  very  small  or  very  large 
vector  elements  is  avoided.  For  example,  the  largest 
vector  element  or  the  vector  length  can  be  normalized  to 
unity  (  5 : 145 ;  7 :  5 . 4). 

The  choice  of  the  vector  v  affects  the  quotient  con¬ 
vergence  rate.  The  quickest  convergence  is  realized  when 
v  is  chosen  such  that 

[x‘*">]Tx“"'>/([x‘k">]TxM) 

and 

the  first  being  the  Rayleigh  quotient  ( 5 : 144 ; 10 : 460 ; 

7:5. 1-5. 8) . 

Kreyszig  (10:462-463)  gives  a  method  for  determining  an 
error  limit  for  the  approximation  Aq  —  Q,  ,  when 

v=x'*  .  Y/here  A«j,  is  the  eigenvalue  closest  to  the 
quotient  approximation.  Note  that  A<^  may  not  equal 
X,  for  early  iterations,  but  will  eventually.  The 

following  is  a  generalization  of  that  method  applicable  to 
both  Aj  —  and  A^  —  and  for  any  arbitrary 

vector  v.  Derivations  are  given  in  Appendix  A.  For 
X^  —  the  error  is  limited  by 

€,  -  S(<i'  -  *■  ^  A7?*  )  2  S, 
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where 


^  /  =  I  ^  1 1 


III 

\/ 

T 

(k+i) 

X 

/^T 

m0  = 

[ 

X 

(k)  jT 

(k) 

X 

r 

('k)'l  T 

(’k  +  0 

P9(  = 

L 

X 

J 

X 

^  s 

[ 

X 

r  (<<♦>) 
X 

For  A,  - 


* 


the  error  is  limited  by 


where 


6£  4  /(cf_l  -  Z  <•  /m.)  =  S£ 

-  I  X%  "  <}.«  I 

*.  «  /TxM/(v'xW) 

w».  «  [  X  <k>]  T  X  (K) 

W,i«txfl')JT  Xft*° 

r  0«*0~1T  (Kti) 

r77+  :  l  X  J  J  X 

For  Xt-  v/^J  ,  the  error  is  limited  by  either 

e*  -  -  /^i)  *■  +  s*) 


if 


or 


£j  -  J(CL^)  ~  J(\i  +  $*) 
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if  X 


y  ~  and  ~  $>1 


,  or 


if  ~  ^ i  and  %*■  ~  *  v,here 

6$  ~  I  ~ 

Iterative  Hatric  Convergence  (10:236-237) 

Recall  that  the  implicit  matric  equation 

> 

A  x  =  b 

was  rewritten  as 

(1)  x  =  G  x  +  c 

and  an  iterative  matric  equation 

(n*0  s~  (n) 

(2)  X  =  G  xl  +c 

was  written  for  the  Jacobi,  Gauss-Seidel ,  and  SOR  methods. 
The  error  associated  with  the  kth  iterated  vector  is  defined 


S  X  -  X(0) 


The  iterated  vectors  will  eventually  converge  to  the 
solution,  x,  only  if  the  error  vectors  converge  to  zero  as  k 
increases  toward  infinity.  A  relation  between  successive 
iterative  errors  can  be  found  by  subtracting  Equation  (2) 
from  Equation  (1)  and  substituting  Equation  (3)  into  the 


result. 


~  r  (h) 

x  -  x  *  G  X  -  (7  X 
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-a'**0  =  6  (  x  -  x  <k>) 


e<kr,)  =  6  6  (k) 


Me  now  consider  those  iterative  matrices  which  have 
eigenvectors  that  span  the  iterative  vector  space.  Those 
include  all  real  symmetric  and  some  real  nonsymmetric 
matrices.  The  convergence  condition  is  found  by  expanding 
the  error  vectors  in  terms  of  the  eigenvectors. 

I  -  I 

6<w=  Gk  eM 

|  (c:c;) 

(*)  e(k)-  l  (c;Gkc,) 

(7)  e(h)  =  2  (c;  A;ke;) 

(h) 

In  order  for  £  to  converge  to  zero  for  all  C ;  and 
C;  ,  the  following  must  be  true: 

lisri  ■>  H 

k-»  oo  A  ;  =  O 


A,  <  I 


for  all  i 


This  is  the  convergence  condition  for  iterative  matrices 
that  have  eigenvectors  which  span  the  iterative  vector  space, 
Unfortunately,  some  iterative  matrices  are  nonsymmetric 
and  have  eigenvalues  of  multiplicity  greater  than  unity.  In 
these  cases  the  eigenvectors  might  need  to  be  supplemented 
by  principal  vectors  in  order  to  span  the  vector  space 
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(5 : 28-30;8 :8 , 17) .  In  these  eases,  some  of  the  terms  in 
Equations  (6)  and  (7)  like 

(9)  C;  Ok  <t;  ~  C;  \  f 


are  replaced  by  terms,  derived  in  Appendix  B,  like 

k-j 


(10) 


C,-6kt:  «  C,  C  t;.-] 


where 


^j)  =  binomial  coefficient  =  k  !  /C\  !  L  U -j  ]  f) 

(3:624) 


=  =  an  eigenvector  associated  with  the 

multiple  eigenvalue 

t;  =  principal  vectors  associated  with 
k  -  i -  rr>  ,  otherwise  for  k  -  )  -  ?Y)  ,  j  =  <7,  k 


It  is  obvious  that  p  <  I  is  a  necessary  condition  for 
convergence  but  it  may  not  be  a  sufficient  condition.  For 
large  multiplicities  and  large  eigenvalues,  some  of  the 
products 


may  not  decrease  to  zero  as  k  increases  to  infinity. 

Mote  Equation  (9)  is  just  a  special  case  of  the  more 
general  Equation  (10)  with  j=0  signifying  that  no  principal 
vectors  are  needed  to  span  the  vector  space. 

If  the  kth  iterated  error  vector  is  examined 


(k) 

e 


I  ( q  a;  e,) 
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or  more  generally 


«"  •  i  U,  %  [  (!)  x.~'  t,.j ) ) 


the  size  of  the  eigenvalues  can  be  seen  to  effect  the 
convergence  rate.  Figure  1  shov/s  that  as  k  increases,  X** 
decreases  faster  for  smaller  X  .  The  spectral  radius, 
being  the  largest  eigenvalue  in  magnitude,  gives  the  slowest 

and  is  therefore  used  as  an  index  for 

convergence  rate. 

Care  must  be  taken  when  comparing  matric  convergence 
rates  by  comparing  spectral  radii.  It  is  quite  possible 
for  the  matrix  v/ith  the  larger  spectral  radius  initially  to 
converge  faster.  The  matrix  v/ith  the  smaller  spectral 
radius  may  have  much  larger  secondary  eigenvalues  than  those 
of  the  matrix  with  the  larger  spectral  radius. 
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Also  the  choice  of  initial  vector  may  be  linearly  inde¬ 
pendent  of  the  eigenvector  associated  v/ith  the  spectral 
radius . 
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> 
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Comparing  two  matrices  v/ith  the  same  spectral  radius,  the 
one  v/ith  the  largest  secondary  eigenvalue  will  eventually 
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PO'.rers  of  Li.r'enval ue 


converge  slower. 


[o  :][:]• 


The  following  discussion  is  limited  to  the  matrix 
problem 

A  x  =  b 


where  A  is  real,  symmetric,  and  block  tridiagonal  (5:45-46). 
The  discussion  can  be  generalized  for  matrices  with  other 
consistent  orderings  (8). 

The  Jacobi  iteration  matrix  is  real  and  symmetric  and 
so  its  eigenvalues  are  real.  Due  to  the  consistent  ordering 
of  the  elements  of  matrix  A,  the  nonzero  eigenvalues  of  the 
Jacobi  matrix  form  positive  and  negative  pairs  (8).  Also 
due  to  the  consistent  ordering  of  the  elements  in  matrix  A, 
a  special  relationship  exists  between  the  eigenvalues  of  the 
Jacobi  and  SOR  iteration  matrices.  (Figure  2) 

(A  ♦  W  -  0*  =  A  (12:243) 

X  ♦  W~  I  =  /(X)  W  jJ 

\(0  +  A*  A vvyw)  *■  (vv-l)  =  o 

XVl  =  (w  fA  t  /Cwy-  ±(\N-  l)J  )  /  Z 

\  -  (wfj/2.  -  /f w*yu,-f(w-03/2)  (8) 

- ( wfj / z 
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bi  Eigenvalues  and  Acceleration  Constant 


X  “  (  w  fA  /  2  -  /  C  (  WyU  /Z^-W+ij) 

=  (  WyH  ll  )  2  +  l(  w /W  /2  )  *  *  /  -  W  3 

i  w p  J L  (  WfA  /2.)Z  +  I  -  wj 
=  w  Ifj1  /  Z  +  I  -  w  ±  WyU  J  ( /  $  *  I  -  W  ) 

where 

yK  =  Jacobi  eigenvalue 

vV  =  SOR  acceleration  constant 

A  =  SOR  eigenvalue 

I7e  now  note  some  basic  properties  of  this  relationship 

(3). 

1 .  No  A  —  O 

(  A  +  w  ~  i  )*  ~  Aw  * 

A  =  o 

(  \A/  -  i  )  2  =  0 

vv  -  |  =  o 

w  =  I 

but  w  >  I  for  SOR. 

2.  The  only  negative  real  eigenvalue  is  A  =  I  ~  W 

3.  X  «=  I  -  vv  if  and  only  if  p  -  O 

2b 


A  *  w  -  I  -  o 
X  -  I  -  w 

and 


(A  t  w  -  l)  =  A  W 1 1*1 

X  =  /  -  w 

0  -  ( I  -  w) 

fJ  *  o 

O  -  ( I  -  w)  vv  1 

W  —  I  or  W  =  O 

but  w  >  /  for  SOR 

4.  A  is  positive  and  real  if  the  radical  is  greater  than 

zero.  I  pi  is  then  >  2,  J(  w-i)/vV. 

W  1  / A'  *  I  -  VV  >  O 

w  1/4 2  /  4-  >  VV  “  / 

/v  1  >  4- 
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\^\  >  2  y<fw-/)/w 

5.  X  =  v/-l  if  and  only  if  the  radical  equals  zero 

X  =  W* fAX / Z  +  I  “  W  -  Wf4  J (  vv  1  /  4-  *  I  -  w) 
w*" yiz  /  +  +  \  -  w  -  O 

-  \N*~  /  1  +1  -W 

W  lJL4Z  /  4-  +  I  ”  W  -  O 
Wl^f  V+  -  w  -  / 
w 1  /  Z  -  2.  \a/  -  Z 

-  2  w  -  Z  +  /  -  vv 

-  w  -  / 

and  I  I  =  Z  J  ( iv  -  I  )  /  vV 

VV  ^  /  4-  +  /  -  vv  =  (9 

*/4-  -  w  -  / 

-  4-  w  -  i)  /  W 
\}aI  -  z  J(w-i)  /  w 

6.  The  number  of  positive  eigenvalues  >  W-/  =  the 

number  of  positive  eigenvalues  <  vv-/  . 

7.  X  is  complex  if  and  only  if  the  radical  is  less  than 

zero  and  £  O  •  If* I  is  then  <  2/^  vv  -  i)/  vV. 
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t  I  ~  W  <  o 

w‘^74  <  vv  -  / 

^l<  4  (w-0  /  vv  * 

1^/  *  2  J(w-i)  /  w 

8.  If  X  is  complex,  the  modulus  of  A  equals  v,--l  . 

J(X)  =  W  fj/ 2.  t  /[  w  /y1*  -  4  ^w-  /)]  /  Z 
W  2yW  ^  -  4-  (  VV _  l)  <  O 
=  WyW/2  f  )  y [4-(w-0  -  WZ{Al  J  / Z 
IX I  =  J(  A*A)  *  [/(A)]*  /(A) 

=  wz//t/f2/4-  +  £  4  ( w-  0'w1’(^iJ/4- 
=  W*~  fA1,  /  4-  +■  vv  -  I  -  vv  *  JAU  /  4- 
-  vv  -  / 

The  mapping  betv/een  Jacobi  and  SOR  eigenvalues 

A  =  w  /*z  /  Z  +  \  -  w  i  vs  fA  J (ssl fA  V4  +  I  ~  vv,) 

gives  a  pair  of  unequal  SOR  eigenvalues  for  each  pair  of 
Jacobi  eigenvalues  unless  the  radical  equals  zero  or  unless 
ra=0.  When  the  radical  equals  zero  and  n^O,  then  the  Jacobi 
eigenvalue  pair  naps  into  a  pair  of  equal  SOR  eigenvalues. 

A  Jacobi  eigenvalue  equal  to  zero  gives  a  one  to  one 
napping. 

A  special  relationship  betv/een  eigenvalues  of  the 


2S 


Jacobi  and  Gauss-Seidel  iteration  matrices  is  obtained  as 
a  limiting  case  of  the  Jacobi-SOR  relationship  where  u=l 

X  -  w Xmfjt / Z  *  I  -  w  i  w J(  w ^  /4  *  I  -  w) 

w  =  | 

(li)  A  *  p  x/&  i  p'/z 

Note  that 

1.  X  is  always  real 

2.  half  of  the  A '-s  -  O 

3.  the  other  half  of  the  Ai  *  O 

Note  that  the  minimum  SOR  spectral  radius  with  respect  to 
the  acceleration  constant  occurs  when  the  radical  equals 
zero,  therefore 

W  Zy£/f  2  /  4-  +  /  -  W  -  O 

vs/* [a1  / A-  -  w  -  I 

fyi  2  =  4  (w-i)  /  vv  * 

=  4  /w  -  4  /  w2" 

/  —  yU  1  =  f/w1  -4/w+| 

=  (e/w  -  <)* 

/^/  -/ut)  =  z/w  -  i 
I  +  J(\ -Hl)  =  z/w' 
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w  =  i  /  L  i  +  J  ( /  ~ 


wb  =  2  /  tit  VO  -  jo(&)  )]  (8  ;  12: 2  A3 ) 


where 


r 


(e) 


=  Jacobi  spectral  radius 


Wb  =  optimum  acceleration  constant 


Thus,  an  optimum  acceleration  constant  (Figure  3)  exists 
that  minimizes  the  spectral  radius  and  so  maximizes  conver¬ 
gence. 

The  spectral  radius  corresponding  to  the  optimum 
acceleration  constant  can  be  found  by  recalling  that  v/hen 
the  radical  equals  zero,  A  =  w-1. 

p(sofC)  =  wk  -  I 

Wt=  2/(1  *J[ i-f(eri)- 1 

(i2)  =2/(1*  /[ pc&yD- 1 

Figure  4  compares  the  Jacobi,  Gauss-Seidel ,  and  optimum  SOR 
spectral  radii. 

Again,  due  to  the  consistent  ordering  of  the  elements 
in  matrix  A,  a  special  relationship  exists  between  the  eigen¬ 
vectors  of  the  Jacobi  and  SOR  iteration  matrices  (8). 

y  -  E  L /( A) ]  V 


where 


y  = 


SOR  eigenvector 


A 


SOR  eigenvalue 
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•;ure  3.  Optimum  Acceleration  Constant 


GAUSS- SEIDEL 


Figure  4.  Spectral  Uadius  Comparison 


m 


V  =  Jacobi  eigenvector  associated  with  the  Jacobi 
eigenvalue  -  O  that  maps  into  X 


E[*J 


a  diagonal  matrix  in  which  the  identity  sub¬ 
matrices  Ij  are  of  the  same  order  as  the  diagonal 
submatrices  of  the  block  tridiagonal  matrix  A 
(5:46)  . 


If  the  radical  ^  0,  then  the  Jacobi  eigenvectors  will 
give  an  equal  number  of  linearly  independent  SOR  eigen¬ 
vectors.  But  if  the  radical  =  0,  then  the  Jacobi  eigen¬ 
vectors  will  yield  only  half  as  many  linearly  independent 
SOR  eigenvectors.  Principal  vectors  will  have  to  be  added 
in  order  to  span  the  vector  space  (8). 

p  =  f  [  y  ( tv  -  i)j  v  /  y c  w~  * 3 

where 

p  =  principal  vector 

V  =  Jacobi  eigenvector  associated  with  the 

that  maps  into  A  s  W -  f 

F  L*J  =  the  diagonal  matrix  with  elements 

f*j  =  cl  E  j  /J  x 

The  special  relationship  between  eigenvectors  of  the 
Jacobi  and  Gauss-Seidel  iteration  matrices  is  obtained  as  a 
limiting  case  of  the  Jacobi-SOR  relation  where  w=l.  The 
relation  is  unchanged  (8) 


34 


E  [  JCx) ]  V 


y  = 

where 

y  =  Gauss-Seidel  eigenvector 
\  =  Gauss-Seidel  eigenvalue 

v  =  Jacobi  eigenvector  associated  with  the  Jacobi 
eigenvalue  fA  -  O  that  maps  into  A 

E  L  =  Same  as  previously  described 

Mote  that  Gauss-Seidel  eigenvectors  are  real  and  span  the 
vector  space. 

Generally  the  Gauss-Seidel  and  SOR  matrices  are  not 
symmetric  so  that  these  special  relationships  are  welcome. 

norms  (6:171-173,176) 

It  is  often  convenient  to  measure  the  size  of  a 
vector  or  a  square  matrix  by  assigning  a  scalar  number  to 
it.  A  norm  does  just  that.  This  section  introduces  the 
concept  of  norm  which  will  be  used  later  to  measure 
iterative  vector  ei’rors,  The  norm  will  be  defined  and  some 
useful  relationships  given. 

A  norm  is  symbolized  by  placing  W  II  around  the 
vector  or  matrix  symbol.  Properties  that  define  a  norm 
are  analogous  to  properties  that  define  an  absolute  value. 
In  fact,  an  absolute  value  is  a  norm. 

II A II  *  o 

II  All  =  o  if  and  only  if  A=0 

II  «  A II  =  H  II AH 
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11  A*  B  II  ±  II A  II  +  II  B  II 

A  matrix  norm  must  also  have  the  following  property 

II A  BU  *  II A II  II  an 

A  vector  n-norm  can  be  defined  as 

tvlln  «  7(  f  I  ViD 

where 

V  =  vector 

V . =  ith  vector  element 
1 

Commonly  used  vector  norms  are  the  one-norm,  the  two-norm 
or  euclidean-norm,  and  the  infinite-norm,  respectively 

II  vll,  •-  I  I  Vi  I 

I 

u v«*  -  J(llv-.l') 
i  vi.  -  T  K-l 

Vector  norms  can  be  ordered  as  follows 

Ivt,  i  llv/lj  *  llvll,  ■■■  *  II  v //„ 

A  matric  n-norm  can  be  defined  as 

II  Mll„  =  T'  ( ZlM  vl„  /  II  vBn  ) 

Matric  norms  are  generally  difficult  to  calculate.  Com¬ 
monly  used  matric  norms  are  the  one-norm  and  the  infinite 


36 


norm. 

ii mii,  ■  "T  (z 

HMD-  -  T“  Ul".jl) 

J  J 

where 

HO  , 'j  =  the  rnatric  element  of  the  ith  row  and 

jth  column. 

If  a  matrix  is  symmetric,  then 

II  Mil,  -  II  Mff„ 

The  spectral  radius  is  less  than  or  equal  to  all  rnatric 
n-norms  ( 12: 237) . 

p  -  II  Mlln  for  all  n 

Conte  and  DeBoor  state  without  proof  (6:228),  that  for  any 
£>0  a  rnatric  norm  exists  such  that 

II  Mil  *  +  £ 

The  triangle  inequality  property  of  a  norm  can  be 
extended  as  follows 

UAH-ten  I  *  IIAII-IIBII 

One  last  inequality  can  be  written 

(Ml  11*11  *  II M  v  II  ±  llvll /ll  M-’ll 
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Iterative  Error  Limits 

The  concept  of  contractive  matrices  is  introduced  to 
aid  in  measuring  the  error  between  the  metric  solution  and 
the  iterated  solution.  The  contractive  constant  and  dis¬ 
placement  vector  will  be  introduced  to  estimate  the  error 
vector. 

Recall  the  implicit  metric  equation 

A  x  =  b 


which  v/as  rewritten  as  (12:233) 

x  =  G  x  ■*  o 

so  that  an  iterated  solution  could  be  attempted  using 

x(n,0  ,  q  f  c 


where 

X  =  rnatric  solution 

(o) 

X  =  nth  iterated  solution 

G  =  iteration  matrix 
C  =  known  constant  vector 


The  error  vector  v/as  defined  as 


ew  ;  x  -  x<n) 


and  the  displacement  vector  V/&S  defined  as 

d(n)  Z  x(^°  -  x  Cn) 
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=  G  (V"-0  -  tCn}) 
*  G  e  <''”0  -  6cw 
=  G  e'-0 

(»>)  fo) 

=  e  -(ye 
=  Cl-G) 

O)  een>  =  (1-0 J™ 


This  might  suggest  that  by  knowing  the  displacement  we 
might  be  able  to  find  the  error,  or  at  least  an  approx¬ 
imation  or  upper  limit. 

If  our  iterative  matric  function 

2  (z)  ~  G  2  +  c 

is  contractive,  then  for  all  z,  a  contractive  constant, 
K,  exists  such  that 

O  ±  K  <  I 


and 


II  g  (*)  -  g  (y)U  ~  K  II  x  -  y  || 

for  some  matric  norm  (6:223).  This  defines  a  contractive 
function.  Simple  relationships  between  the  displacement 
and  error  vectors  are  found  to  exist  by  using  this  con- 


(8) 
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tractive  prooerty. 


x  ~  (3  x  +  <s 
x  =  6  x'"'  »  c 
lx-  xf~,>|  6  K  Ix-x^’l 
e4’*  5  x  -  xf,,> 

(4)  II  e"”"|  A  K  Ic'-’l 

x"”1’  = 

X  -ox  +6 

|U <»,o_  A  *”*>|  x  Kllx"”0  -x^’ll 

I  £»}  <*♦»)  O'*) 

d  5  X  -  X 

<5)  m‘-')ii  *  k  kjmh 

x  -  <3  x  +  C 

X  *  Sxw»c 

lx  -  xw>J  K  i  II  X-  X^l 

=  //x-x'”>  ♦  X -  X 

±  II  X  -  X  Mll  -  II  x'"'0  -  x^ll 
lx""’  *  Ix-x'”’!  -I*-X'”IK 
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II  x  -  x II  -  llx-x,n>s  (i-k) 
jM  .  xr». o  .  x  CJ 

eM*  x-  xfo) 

U°"t  i  lie.™  II  (i-  k) 

I  >  K 

I  -  K  >  a 

lleMH  i  lid ‘-’ll  /0-k) 


These  are  analogous  to  the  previous  relationships. 

A  variety  of  relationships  between  errors  and  dis¬ 
placements  of  differing  iterations  can  be  obtained  by 
applying  either  or  both  Equations  (4)  and  (5)  to  (G) 
successively.  For  example, 

(7)  K/(i-k)  (6:223) 

(8)  lld^lt  Kn/(l-K)  (5:223) 

The  same  can  be  done  with  Equations  (1),  (2),  and  (3) 

ew  =  (1-6)"  G  2'-° 
eM  -  (r-6)"  G”  2'° 


Ilote  that 

hMs 4  *  II  /""’l/ «/(/-«)  *  W"’i if  H-yci-w) 

How  is  it  knovm  that  a  function  g(s)  is  contractive? 
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!Io\;  is  it  known  that  a  contractive  constant  less  than  unity 
exists?  Assume  the  rnatric  function  is  contractive. 


II  36*)  -  3  (y)ll  ~  Kllx-yll 

3  =  G  x  +  c 

d(y  )  =  6  y  +  c 

~3^)0-Gx-Gy 
=  G  (x- y) 

V  =  X  -  y 

=  G  v 


II  G  v  II  -  K  II  V H  for  all  v 

II  G  v  II  /  II  v/j  ^  K  <  I  for  all  V 

"v*  II  Gv  II  /II  VII  6  K  <  I 

II  G II  S  II  Ov  ll/llvll 

II  Gil  *  K  <  / 

Ilote  that  II G II  is  always  an  acceptable  A  if  l/G/l 
is  loss  than  unity. 


3  (*■)  -  3^ y)  =  G  (  x  -  y) 

II  $(*)  -  g(y)  H  =  II  6  (x  -  y)ll 

*  II  Gil  II  x  -  y  I/ 


Thus  our  iterative  matric  function  is  contractive  if  and 
only  if  some  matric  norm  of  C.  is  less  than  unity  (6:225). 


116//  <  / 


Recall,  from  the  section  on  norms,  that  a  matric  norm 
exists  such  that 

(g)  -  II G  ll  -  (&)  +  £ 


v;here 


(G)  =  spectral  radius  of  the  matrix  G 

£  =  any  number  greater  than  zero 

Thus  a  matric  function  is  contractive  if  and  only  if  the 
spectral  radius  is  less  than  unity  (6:228) 

<  I 

This  is  the  same  condition  for  convergence. 

Normally  one-norms  and  infinite-norms  are  used  for 
analysis  because  they  are  easy  to  compute  and  they  best 
represent  the  type  data  to  be  analyzed.  Unfortunately, 
the  contractive  relationships  are  valid  only  for  those 
vector  norms  associated  with  matric  norms  less  than  unity. 
Recall,  from  the  section  on  norms,  that  if  the  magnitude 
of  any  matric  element  is  greater  than  unity,  then  the 


inf inite- 

norm  and 

one 

-norm  are  also  greater  than  unity. 

Even 

though 

the 

relationships 

(4) 

1  1 

K  (fee’ll 

(5) 

II 

I  * 

K  lid*”  It 

(C) 


Ve^ll  *  ltdl"",ll  K/O-X) 


(7) 

may  be  invalid  for  K  limited  to  1  >  K  ±  0,  extended 

domains  of  K  exist  where  the  relationships  are  always  valid, 
whether  the  matric  function  is  contractive  or  not. 

Recall  Equation  (1) 

eC""')  =  Q  d« 

II  e |  =  |  6c  c”'l 

^  II  G II  lie"" II 

Comparing  this  with  Equation  (4)  shows  that 

O  ±  K  ±  II G  II 

is  the  extended  domain  of  K  for  Equation  (4).  Similarly, 
recall  Equation  (2) 

=  6  /•> 

II  J(”">ll  *  II  G  JMII 

*  Hell  lids’ll 

Comparing  this  wi th  Equation  (5)  shows  that 

O  -  K  ^  I/O  II 

is  also  the  extended  domain  of  VC  for  Equation  (5).  By 
inspection,  the  extended  domain  of  IC  for  Equation  (6)  is 
(Figure  5) 
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\ 


-<*>  6  ft  <  I 

By  inspection,  the  domain  of  K  for  Equation  (7)  needs  no 
extension  (Figure  5) 

O  *  ft  <  I 


Finally  v/e  analyze  hov/  these  four  relationships 
behave  in  the  limit  as  the  number  of  iterations  increases 
toward  infinity.  This  is  the  region  where  our  iterative 
solution  converges  to  our  matric  solution.  Note  that  the 
following  derivations  assume  that  the  eigenvectors  span 
the  vector  space  and  that  the  principal  eigenvalue  is 
real  and  does  not  have  a  negative  counterpart 


A,  ¥■  -  A; 


for  any  i 


Recall  the  sections  on  spectral  radius  and  iterative  matric 
convergence  for  difficulties  that  can  occur. 

First  we  derive  the  j'^oo  lle(ntlill  /  lle.(n>ll 
Expanding  the  initial  error  in  terms  of  the  eigenvectors, 
v,  of  the  matrix,  G  (12:238) 

=  Z  (c,v.) 

l-l  ' 

(n)  r  n  (o) 

e  -  C?  e 
=  G"  I  (C;V;) 

l-l 

=  |  (c;  6"  Vi) 

=  fc.U/AjV,)] 
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is  of  multiplicity  s 


’./here  X 


X;  <  A, 

I  A;  /A,  I  <  I 

[A; /A  ,r=o 
w  =  a;  i(cjv;) 


III*) 
r\  -*  eo 


I,' i*) 

n  -*oo  S. 


v/hich  is  an  eigenvector. 

(10) 


/.V77  (o  n)  _  v 

n-*o°  S  "  -»«o  £  A  t 


Using  the  quotient  method  described  in  the  section  on 
spectral  radius 

T  (<Otl)  /  /  T  Co)  \  _  V 

we.  /  (  vV  e  /  -  A, 

or  using  norms 


HZ-  Ic^0!  =  I  A,  e-^fl 
=  /  A,/  'tu  He"”// 


(11) 

birj 

II  e 

’ll/ Ilex’ll  -  IX,  1 

III 

(&) 

Similarly 

(10: ; 

41) 

hrn 
n-t  *0 

1/rO 

=  n  -*e° 

J(n)  A, 

hnO 

0  -+CO 

W  T  c 

\in"’  /  ( 

wTo/^)  «  A, 
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(12)  /UMH  «  /  A,/ 

Next  we  derive  n-r*»  He.*  *11  /lid*  II 


(13) 


(n)  (*) 

e  -  X  -  X 


(an) 


(an) 


(*) 


e  '  *  x  -  x 

cw.  t <•»..) ,  *<»♦■>  -  x 

d"0--  x^1 

e  -  c 

)»'*7  /  (a)  (nil)  \  _  /»>>)  l(n) 

n  V.  C  ~  C.  J  ~  n  O 

/.'•»>  /  fo)  »  (a)\  _  />*vj  j(n) 

n~toQ  \  £■  —  Aj  £  /  ~  0-9*0  0 

I . 


liar, 

n  -*  co 
h*n 


S"(<  -  A,)  •  «/ 


nZ«  lle‘n'0-*,)h 

ll  *‘"11  1 1 -a,  I  = IIS"  II 

IIS"!/ IIS"  II  =  I //I-A,  I 


I  <*n 
n  -9bo 


ll’l'n 

o 


x,  <  / 

IIS"  11/ US" II  =  i  /  (  I  -  x.) 


Mote,  the  difficulty  of  having  A,  *  -  A;  is 
from  this  last  equation.  Two  limits  would  exist, 
nonsense . 

Similarly 
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( 8 ; 1 2 : 240) 


apparent 
which  is 


.  —  i-V 


/.>»>  (n)  /  .  \  \  _  i  (n  - 1 )  . 

n-**°  e  (.  I  ~  A, )  -  a  A , 

(14)  i*.  lle.^11  /II  =  /X,/  /(>- A,) 

Note  that  X,  <  o  can  not  be  a  K  in  the  limit  because  of 
the  absolute  value.  An  equivalent  K  can  be  derived. 

(Figure  6) 

"  k/6-k)  =  IX,I  /  ( (  -  Aj 

K  (t- A,)  =  IX.  I  (i-k) 

K  -  K  A,  -  /  A,  |  -  K  |  A,  I 
K  -  K  A,  *  K  l\,l  =  IX,  | 

K  0-X,  +  I  A,  I )  =  IX, I 
(is)  K '  IX.l  /  (  I-  X,  +  t\.l) 

If  the  principal  eigenvalue  is  negative 

x,  <  o 

A.  =  -  /  A,  I 

K  *  IX,  |  /0  ♦  I  A,  I  -  I  X.l) 

(16)  -  I  A,  |  /  (l*  2IX.I) 

If  the  principal  eigenvalue  is  positive 

A,  >  a 

A,  =  lx,  I 
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Figure  6.  The  K  that  satisfies  h'n  He<'>*ll/ lld(n~'*l=  K/(I-K) 


K  =  IX, I  /  (  I  ~  IX,  I  ‘  IX, 
=  IX,  I 

as  expected. 
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III.  Procedure 


Sample  Problem 

Me  now  pose  the  sample  problem  that  will  be  the 
subject  of  experimentation.  A  parabolic  two  dimensional 
partial  differential  equation  was  written  to  represent 
transient  heat  conduction 

(i/k)  >T/Jt  y t/J*'  <-  yr/iy* 


It  was  then  approximated  using  a  full  implicit  form  of 
differencing 

( |/K )  A,  T,- „  / 

=  CT;^,  /C  *  j  /h/ 

’which  was  then  expanded 

(t:-  „  -  T,  „ )  /(hf  K ) 

+  .  -  ZT.  .  <-  T  V 

terms  collected  with  unknowns  on  left  side 

Tw>  =TW,~.  0*2<,,kA/  .Zh.K/h,‘) 

«■  (-KK/hS)*  (- h, kA;) 

♦W<  K*/V> 


and  constants  defined 
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+ 

b 

T-  • 

+• 

c 

T.  • 

V  J  /'■>#■/ 

+ 

b 

T-  .  „ 

»/  j  +  i,  "'i 

+■ 

a 

T.  .  . 

where 

a  =  -Kb,/  h* 
b  =  -  K  lof  / hyl 
c  =  I  *  2  hf  K/ hxl  +  2  ht  K / 

The  computational  molecule  (12:12)  that  represents 
the  relationships  between  variables  of  the  difference 
equation  is  depicted  in  Figure  7.  The  special  subscripts 
must  be  combined  before  the  function  T  can  be  represented  by 
a  vector.  The  interior  nodes  are  equated  to  sequential 
vector  elements.  If  the  spacial  domain  is  rectangular, 
the  subscript  i  is  incremented  through  its  entire  range 
while  for  each  increment  of  i,  the  subscript  j  is  incre¬ 
mented  through  its  entire  range.  The  resulting  computational 
molecule  is  depicted  in  Figure  8. 

This  sequencing  of  interior  nodes  results  in  an 
implicit  matric  equation 

Ax  =  b 

with  the  matrix,  A,  being  block  tridiagonal.  The  sub¬ 
matrices  on  the  diagonal  will  be  tridiagonal  with  diagonal 
elements  equal  to  c  and  off  diagonal  elements  equal  to  b. 

The  submatrices  on  the  off  diagonal  will  be  diagonal  with 
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Figure  7.  Computational  '.iolecule  Usinj  T’.;o  Spacial 
Indices 


n+1 

n 


Figure  3.  Computational  liolecule  Usinp  One  Special 
Index 


diagonal  eleiaents  equal  to  a.  For  example,  a  rectangular 
spacial  domain  with  two  interior  nodes  on  a  side  will 
result  in 

C  b  q  o 

be  o  a 
a  0  c  b 
0  a  b  c 

Appendix  C  contains  the  analytic  solution  to  the  4  by  4 
implicit  rnatric  equation,  contains  the  Jacobi,  Gauss-Seidel , 
and  SOR  4  by  4  iterative  rnatric  equations  ctnd  the  associated 
eigenvalues.  This  information  was  for  computer  code 
validation. 

It  was  necessary  to  pick  initial  and  boundary 
conditions  that  resulted  in  a  partial  differential  equation 
that  could  be  solved  analytically.  This  physical  solution 
could  later  be  comoared  to  rnatric  solutions,  x,  and  to 

/  n ) 

iterated  solutions,  x'' 

The  analytic  solution  of  the  parabolic  partial 
differential  equation 

J'T/Jx’  *  i'r/i  yL  =  (</ k)  J T/it 

over  the  spacial  and  time  domains 

o  ±  x  ±  ir 

0  £  y  6TT 
O  6  t 

with  homogenious  Dirichlet  boundary  conditions  (1:333) 

T (o,  y,  r)  =  T(ir,  y,  t)  -  o 


T (*/o,is)  -  T  ( xy  TT/  t)  -  O 

and  an  initial  condition 

T  (*,  y  ,  o  )  x  5  i  *o  (.*  )  t>m  (y) 


is  equal  to 


-2  KT 

T(\y,t)--  Sin  (x)  s in  ( y >  « 


The  derivation  and  validation  is  contained  in  Appendix  D. 

Being  that  the  boundary  conditions  are  homogeneous, 
the  vector  b  of  the  implicit  metric  equation 


Ax  -  b 


is  equal  to  the  initial  conditions.  The  homogeneous 
boundary  conditions  contribute  only  terms  equal  to  zero. 

Computer  Codes 

Two  computer  codes  were  designed  to  measure  the 
iterative  metric  spectral  radius  as  a  function  of  the 
density  of  nodes  over  the  special  variable  domain.  The 
Jacobi  spectral  radius  v/as  first  determined  using  the  power 
method,  then  the  Gauss-Seidel  and  SOR  spectral  radii  were 
calculated  from  Equc-tions  (11)  and  (12).  The  first  code 
determines  the  spectral  radius  using  a  quotient  of  iterative 
vectors  separated  by  one  iteration  (Appendix  E) ,  the  second 
uses  a  quotient  of  iterated  vectors  separated  by  two 
iterations  (Appendix  F).  The  second  code  was  run  to  make 
certain  that  truncation  error  forced  the  two  maximum 
modulus  Jacobi  eigenvalues  to  have  slightly  differing 
magitudes. 

One  computer  code  (Appendix  G)  v/as  designed  to  iterate 
sample  problem  solutions  by  using  successive  over-relaxation. 
Analytic  rnatric  solutions  that  corresponded  to  physical 
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solutions  could  only  be  calculated  for  the  low  nodal 
densities,  3  by  3 ,  4  by  4,  and  5  by  5.  liatric  solutions 
for  higher  nodal  densities  were  not  correlated  with  physical 
solutions.  In  these  cases  an  initial  iterative  vector  was 
calculated  that  would  converge  to  a  preplanned  rnatric 
solution  by  solving  the  rnatric  equation  in  reverse. 

A  x  =  b 

where 

X  =  preplanned  matric  solution 

b  =  calculated  initial  iteration  vector 

In  this  way  the  error  between  the  iterated  vector  and  the 
rnatric  solution  was  easy  to  calculate. 

Five  sets  of  data  were  generated  for  later  comparison. 
The  first  set  of  data  contains  the  physical  and  rnatric 
solutions  and,  for  each  iteration,  an  iterated  solution. 

The  errors  between  each  pair  of  solutions  is  also  given. 

The  second  set  contains  displacement  norm  ratios  and  error 
norm  ratios  which  can  be  compared  to  their  limiting  value, 
the  spectral  radius.  The  third  set  contains  the  values  of 
the  parameter  that  is  required  to  make  each  error  approxi¬ 
mation  an  equality. 

lie ‘"’ll  =  IUMII  /(l-K) 

He^ll  o-k)  =  ilJ^i 

(i-k)  *  /lie ‘"’ll 

K  -  I  -  /lie"’’ I I 
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nMi  K/o-tr) 

U‘"'H  (i-n)  =  « 

lice’ll  *  (  IU‘"-°II  <■  lU'"‘/0  K 
K  -  He  II  /  (  II J  '""'ll  *  II  e  ) 

K  -  '/(  II  lie's'll  <-/) 

They  will  reveal  how  the  parameter  varies  from  iteration 
to  iteration  and  what  bounds  exist,  if  any.  The  fourth  set 
of  data  contains  the  error  norm  to  displacement  norm  ratios 
actually  observed  for  two  of  the  error  approximation 
methods. 


IU‘"'ll  =  SJ‘",II  /  O- k) 
lle‘"‘ll  /  II  J ‘"'/I  =  I  /( i-K ) 

IU‘"'l  *  IIJ'"-'H  K/(l-K) 

lle‘"'ll/ll  =  K/(l-K) 

They  will  reveal  how  the  ratio  varies  from  iteration  to 
iteration  and  what  bounds  exist,  if  any.  The  fifth  set  of 
data  contains  approximations  to  the  error  between  iterative 
and  matric  solution  using  the  following  three  error  methods. 

iu‘"'ll  ±  IU‘"'I  / (i-k) 

II  e '"‘11  *  IUln-‘,l  k/P-kJ 
lie ‘"'II  <  lldMll  Kn/(i-K) 
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v/here  the  parameter  IC  which  satisfies  the  inequality  i 
approximated  by  the  following  six  scalar  quantities 


II  Jl" "’ll /tlJ'",ll  *k 

II  J(n)i  /  IIJ‘n  "’ll  -  K 


(Jacobi)  ~  K 
O  ( Gauss-Seidel )  ^  ^ 


IV  Numerical  Results 


Spectral  Radius  Versus  IJodal  Density 

This  next  topic  concerns  spectral  radius  of  the  itera¬ 
tive  matrix  versus  the  spacial  node  density  and  the  size  of 
time  interval  between  solutions  of  the  implicit  metric 
equation. 

Figure  9  depicts  the  Jacobi  spectral  radius  as  a 
function  of  nodal  density  in  a  square  spacial  domain. 

Figure  10  depicts  the  corresponding  SOI’  optimum  spectral 
radius.  Figures  9  and  10  show  that  the  spectral  radius 
increases  with  increasing  nodal  density. 

Tables  II,  III,  IV,  and  V  list  the  Ga.uss-Seidel 
spectral  radii  as  a  function  of  spacial  nodal  density. 

Each  table  represents  a  different  time  interval  between 
solutions  of  the  implicit  metric  equation.  Similarly 
Tables  VI,  VII,  VIII,  and  IX  list  the  SOR  optimum  spectral 
radii.  Figure  11  graphically  depicts  values  from  Tables 
VI,  VII,  VIII,  and  IX  where  x  nodal  density  equals  y  nodal 
density.  The  tables  and  Figure  11  show  that  the  spectral 
radius  increases  with  increasing  size  of  time  interval. 

The  change  in  spectral  radius  'with  respect  to  the  time 
interval  is  largest  when  the  time  interval  is  small  and 
assymptatically  approaches  zero  as  the  time  interval  gets 
large. 

These  increases  in  spectral  radius  cause  the  itera¬ 
tive  matric  equation  to  converge  slower.  See  the  discussion 
on  page  22  and  Figure  1.  Y/hen  one  tries  to  decrease  the 
truncation  error  between  physical  and  matric  solutions  by 
increasing  nodal  density,  error  between  the  matric  and 
iterative  solutions  is  inadvertently  increased  by  slowing 
the  convergence  rate.  An  example  of  this  is  depicted  in 
Figure  12  where  the  convergence  rate  using  3  by  3 ,  4  by  4, 
and  5  by  5  nodal  densities  can  be  compared.  Each  line 
represents  iterative  approximations  to  the  sample  physical 
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crsus  Nod 


Table  II 


h 


Gauss-Scidel  Spectral  Hadius 
Time  Interval  =  1 


it-* ** 

3* 

4* 

5* 

6* 

7* 

8* 

9* 

3 

.154 

.279 

.403 

.510 

.597 

.667 

.721 

4 

.279 

.377 

.471 

.557 

.629 

.688 

.736 

5 

.403 

.471 

.542 

.  608 

.665 

.714 

.754 

6 

.510 

.  557 

.  608 

.657 

.701 

.741 

.  774 

7 

.  597 

.629 

.665 

.  701 

.736 

.767 

.794 

8 

.667 

.688 

.714 

.741 

.767 

.791 

.814 

9 

.721 

.73  6 

.754 

.774 

.794 

.814 

.832 

* Number 

of  node 

intervals  along 

x-axis 

** Number 

of  node 

interva 

Is  along 

y-axis 

Table 

III 

Gauss 

-Seidel  Soectral 

Radius 

Time  Interval  = 

10 

# * 

3  * 

4* 

5* 

6* 

7* 

8* 

9* 

3 

.23  7 

.385 

.514 

.615 

.692 

.750 

.794 

4 

.385 

.487 

.578 

.  656 

.717 

.765 

.804 

5 

.514 

.578 

.642 

.  699 

.747 

.785 

.318 

6 

.615 

.  656 

.  699 

.740 

.776 

.803 

.834 

7 

.692 

.717 

.747 

.  776 

.804 

.828 

.849 

8 

.750 

.766 

.786 

.808 

.828 

.847 

.864 

9 

.794 

.304 

.318 

.834 

.849 

.864 

.870 

*  Number  of  node  intervals  along  x-axis 

**Number  of  node  intervals  along  y-axis 


Table  IV 


Gauss-Seidel  Spectral  Radius 
Time  Interval  =  100 


** 

3* 

4* 

5* 

6* 

7* 

8* 

9* 

3 

.249 

.399 

.527 

.627 

.702 

.759 

.801 

4 

.399 

.500 

.590 

.667 

.727 

.775 

.812 

5 

.527 

.590 

.653 

.709 

.755 

.794 

.825 

6 

.627 

.667 

.709 

.749 

.784 

.815 

.840 

7 

.702 

.727 

.756 

.784 

.011 

.835 

.855 

8 

.759 

.775 

.794 

.815 

.835 

.853 

.869 

9 

.801 

.812 

.825 

.340 

.355 

.869 

.882 

*IJumber 

of  node 

interva 

Is  along 

x-axis 

lumber  of  noie  intervals  along  y-axis 


Table  V 

Gauss-Seidel  Spectral  Radius 
Tine  Interval  =  1000 


** 

3* 

4* 

5* 

6* 

7* 

8* 

9* 

3 

.250 

.400 

.529 

.528 

.704 

.760 

.802 

4 

.400 

.502 

.592 

.  568 

.728 

.  775 

.812 

5 

.529 

.592 

.654 

.710 

.757 

.795 

.826 

6 

.628 

.668 

.710 

.750 

.785 

.815 

H 

CO 

• 

7 

.704 

.728 

.757 

.78  5 

.812 

.335 

.856 

8 

.760 

.775 

.795 

.815 

.  83  5 

.853 

.870 

9 

.802 

.812 

.826 

.841 

.856 

.870 

.883 

*Nur.ibei' 

of  node 

intervals  along 

x-axis 

**JJunber  of  node  intervals  along  y-axis 
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Table  VI 


SOR  Optimum  Spectral  Radius 
Time  Interval  =  1 


** 

3  * 

4* 

5* 

6* 

7 * 

8* 

9* 

3 

.042 

.032 

.128 

.177 

.224 

.268 

.309 

4 

.082 

.118 

.158 

.201 

.243 

.283 

.321 

5 

.128 

.158 

.193 

.230 

.207 

.303 

.337 

6 

.177 

.201 

.230 

.261 

.293 

.325 

.356 

7 

.  224 

.243 

.267 

.293 

.321 

.349 

.376 

3 

.  208 

.283 

.303 

.325 

.349 

.3  73 

.397 

9 

.309 

.321 

.337 

.356 

.376 

.397 

.418 

* Number 

of  node 

interval 

s  along 

x-axis 

**Number 

of  node 

intervals  alone 

y-axis 

Table  VII 

SOR  Optimum  Spectral  Radius 
Time  Interval  =  10 


** 

3* 

4* 

5* 

6* 

7* 

8* 

9* 

3 

.067 

.121 

.178 

.234 

.28  6 

.333 

.375 

4 

.121 

.165 

.212 

.  200 

.300 

.343 

.38  0 

5 

.173 

.212 

.251 

.291 

.330 

.358 

.402 

6 

.234 

.200 

.291 

.325 

.358 

.390 

.421 

7 

.  28  6 

.300 

.330 

.358 

.386 

.414 

.441 

8 

.333 

.348 

.3  58 

.390 

.414 

.438 

.451 

9 

.375 

.  3  S  o 

.402 

.421 

.441 

.461 

.482 

^Number 

of  node 

interva 

Is  along 

x-axis 

**Number  of  node  intervals  along  y-axis 
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Table  VIII 


SOil  Optimum  Spectral  Radius 
Time  Interval  =  100 


*  * 

3* 

4“ 

5* 

6* 

7* 

8* 

9* 

3 

.071 

.126 

.185 

.242 

.294 

.341 

.384 

4 

.126 

.172 

.219 

.268 

.314 

.356 

.395 

5 

.185 

.219 

.259 

.  299 

.338 

.376 

.410 

6 

.  242 

.268 

.299 

.332 

.3  66 

.393 

.429 

7 

.294 

.314 

.338 

.3  66 

.394 

.422 

.449 

o 

O 

.341 

.356 

.3  76 

.398 

.422 

.446 

.469 

9 

.384 

.395 

.410 

.429 

.449 

.469 

.489 

■“Number 

of  node 

intervals  along 

x— axis 

““Number 

of  node 

interv 

als  along 

y-axis 

Table 

T 

S0H  Optimum  Spectral 

Radius 

Time 

Interval  =  : 

L000 

** 

3  “ 

4“ 

5“ 

6* 

7“ 

8“ 

9* 

3 

.072 

.127 

.186 

.243 

.  295 

.342 

.384 

4 

.127 

.172 

.220 

.  259 

.314 

.357 

.395 

5 

.18  6 

.  220 

.250 

.300 

.339 

.377 

.411 

6 

.  243 

.  269 

.300 

.333 

.3  67 

.399 

.430 

7 

.  295 

.314 

.339 

.367 

.395 

.423 

.449 

8 

.342 

.357 

.377 

.399 

.423 

.  445 

.470 

o 

.334 

.395 

.411 

.430 

.449 

.470 

.490 

*Number  of  node  intervals  along  x-axis 
-"'•“Number  of  node  intervals  along  y-axis 
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■ure  11.  Optimum  SOI?  Spectral  Sadi us  versus  Time  Step  Interval 


a 


Convergence  Rate  versus  Modal  Density 


problem  described  earlier. 

Iterative  Error  Approximations 

The  computed  error  data  is  displayed  in  sequences  of 
nine  graphs  each.  Each  sequence  represents  results  using 
a  different  nodal  density  over  the  sane  special  domain  of 
the  sample  problem.  Graphs  representing  a  domain  v/ith  ten 
nodal  intervals  on  a  side  are  included  in  this  section.  The 
data  was  determined  by  using  the  following  selected 
parameters  and  resulting  constants. 

1.  Humber  of  interior  nodes  =  81 

2.  Time  interval  between  metric  solutions  =  10,000 

3.  Diffusion  coefficient  =  1 

4.  Computer  precision  =  27  decimal  digits 

5.  Humber  of  spectral  radius  iterations  =  50 
5.  Jacobi  spectral  radius  =  .05105 

7.  Gauss-Seidel  spectral  radius  =  .90451 
3.  S0R  optimum  spectral  radius  =  .52736 
9.  Humber  of  iterations  =  200 

Graphs  representing  other  nodal  densities  are  included  as 
Appendices  H  and  I. 

The  first  three  graphs  of  a  sequence  depict  relation¬ 
ships  between  error  limit  parameters,  displacement  norms, 
and  error  norms  when  using 

(c>  Bc'-’fl  4  Ill’ll  /(l-K) 

as  the  error  limit  method.  The  second  three  graphs  depict 
relationships  when  using 

(7)  De<"’l  *  K/0-K) 

as  the  error  limit  method.  The  seventh  graph  depicts 
relationships  when  using 

(0)  I  kVO-k) 
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as  the  error  liir.it  method.  Henceforth  these  methods  will 
be  labeled  one,  tv;o,  and  three  respectively.  Note,  here 
all  norms  are  infinite-norms. 

The  first  graph,  Figure  13  for  example,  shows  the 
error  norm  to  displacement  norm  ratio  related  to  the  first 
error  limit  method.  The  second  graph,  Figure  14  for 
example,  shows  the  parameter  K  needed  to  make  Equation  (6) 
an  equality.  Candidate  parameters  may  be  compared  to  this 
minimum  parameter.  The  third  graph,  Figure  15  for  example, 
depicts  the  displacement  norm  and  a  set  of  parallel  curves 
Each  parallel  curve  represents 

Ill’ll  /O-k) 

where  the  parameter  X  is  chosen  to  space  the  curves  an 
order  of  magnitude  apart.  The  error  norm  is  superimposed 
so  that  an  order  of  magnitude  comparison  can  be  made 
between  the  error  nor;.'.,  tne  displacement  norm  and  candidate 
error  limits  based  on  different  parameters. 

Similarly,  the  fourth  graph,  Figure  16  for  example, 
shows  the  error  nor:;  to  displacement  norm  ratio  related  to 
the  second  error  limit  method.  The  fifth  graph,  Figure  17 
for  example,  shows  tne  parameter  X  needed  to  make  Equation 
(7)  an  equality.  Candidate  parameters  again  may  be  com¬ 
pared  to  this  minimum  parameter.  The  sixth  graph.  Figure  13 
for  example,  depicts  the  displacement  norm  and  a  set  of 
parallel  curves.  Each  parallel  curve  remresen-cs 

k/P-k) 

where  the  parameter  X  is  chosen  to  space  the  curv&s  an 
order  of  magnitude  apart,  here  again,  the  error  no ... 
superimposed  so  that  an  order  of  magnitude  com: arisen  e:  . 
be  made  between  the  error  norm,  the  displacement  nor.;,  and 

can.'!  '.at:;  error  co  based  cn  difforono  para..c  cc-rs . 
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The  seventh  graph ,  Figure  19  for  example,  pertains 
to  the  third  error  Unit  Method.  The  ",rnph  depicts  straight 
lines  representing 

IIJ‘1  K"/0-k) 

v/here  the  parameter  K  is  incremented  by  one  tenth  from 
aero  to  unity. 

(20)  i»3  ti/lKVd-K)] 

-  n  log  (ft)  *  log  l(  d  II  ~  log  ( I  ~  kJ 

This  is  in  the  form  of  a  straight  line 

y  -  rri  x  +  b 

v.-here 

y  -  log  [lid  (oill  K  V  (l  -  k)  ] 
x  =  n 

b  =  log  IIJMII  -  Io3  (i-k) 

rr>  *  log  (K) 

The  error  norm  is  superimposed  so  that  an  order  of  magni¬ 
tude  comparison  can  he  made  betv/een  the  error  norm  and 
candidate  error  limits  based  on  different  parameters  d. 

The  eighth  graph,  Figure  20  for  example,  depicts  the 
displacement  norm  ratio.  It  is  a  candidate  error  limit 
parameter,  but  unlike  the  other  candidate  parameters,  it 
is  not  constant.  The  displacement  norm  ratio  can  also  be 
compared  to  its  limiting  value,  the  spectral  radius . 
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1 

O'l 

8*0 

9  '0  t'O  Z'Q 

O’O 

(I+N3/N0)  /I->i 

OUP/Jt? 


of  Magnitude  Comparison  between  Krror  I.'orm  and  Displacement  Norm 
Limit  llethcl  2,  10  ;;  10  Ilodal  Density 


atio,  10  x  10  Modal  Density 


METHOD 


Figure  21.  Three  Error  Limit  Methods  Compared  to  Error  Norm,  K  =  SOU  Optimum 
Spectral  iiadius,  10  x  10  Nodal  Density 


*7V«**Mi*»W*'*>r,  ... 


The  ninth  graph ,  Figure  21  for  example,  compares 
the  error  norin  to  the  three  error  Unit  methods.  The 
parameter  K  in  all  error  limit  methods  being  the  same. 

In  comparing  error  limit  method  one  to  error  limit 
method  two,  it  can  easily  be  shown  that  when  the  dis¬ 
placement  norm  ratio  (Figure  20)  is  less  than  the  parameter 
K 

U(-‘t  /  |N",'°/I  <  K 

then  method  one  is  less  than  method  two 

<  8 k/O-k) 

v.rhile  v;hen 

IIJ^I/l >  k 

method  two  is  less  than  method  one 

I Utn-°l  k/o-k)  <  IIJ^I/O-k) 

Then  the  ratio  equals  ti'.e  parameter, 

«J‘"’l/uin'0S  -  K 

both  methods  give  the  same  result 

UcJ‘'-°ll  k/o-k>  II  J(n>l /(i-k) 

This  does  not  imply  which  method  approximated  the 
error  more  closely.  That  depends  on  which  candidate 
parameter  is  used  in  ea.ch  method.  Figures  14  and  17  show 
the  parameters  that  would  make  each  approximation  an 
equality.  The  following  candidate  parameters  for  the  10 
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by  10  nodal  density  can  be  compared  with  the  graphs  in  this 
section. 

1.  Displacement  norm  ratio:  Figure  in 

2.  Jacobi  spectral  radius  =  .951 

3.  Gauss-Seidel  spectral  radius  =  .905 

4.  SOU  optimum  spectral  radius  =  .528 

5.  Zero 

dote  that  on  all  the  graphs  representing  10  by  10  nodal 
density,  computer  truncation  error  overrides  any  further 
iterative  convergence  past  iteration  number  130. 

Figures  14  and  17  show  that  the  optimum  error  limit 
parameter  oscillates  about  a  value  close  to  the  50R  optimum 
spectral  radius.  The  fact  that  the  SOM  optimum  spectral 
radius  is  the  limiting  value  of  the  optimum  M ,  maxes  it  a 
good  choice.  The  Jacobi  and  C-auss-Seidel  spectral  radii, 
in  this  case,  fall  above  all  optimum  values  of  M.  Using 
a  parameter  equal  to  zero  in  method  one  or  equal  to  one 
half  in  method  two  essentially  approximates  the  error  by 
the  displacement.  It  is  neither  a  lower  limit,  upper  limit, 
nor  as  good  an  approximation  as  when  using  the  30R  optimum 
spectral  radius. 

Figures  15  and  13  show  order  of  magnitude  comparisons 
between  the  error  norm  ana  approximations  using  different 
parameters.  IJote  in  particular  that  as  K  approaches  unity, 
the  order  of  magnitude  changes  rapidly.  These  two  error 
limit  methods  are  very  sensitive  to  changes  in  parameter 
when  the  parameter  approaches  unity,  in  the  limiting  case, 
when  the  spectral  radius  approaches  unity. 

Mote  that  the  error  nor:;*  roughly  parallels  the  dis¬ 
placement  norm,  as  it  should  in  the  limit  as  the  number  of 
iterations  increase.  Doth  lines  have  a  slope  approaching 
the  logarithm  of  the  spectral  radius.  Similar  to  Equations 
(11)  and  (12) 

-  IX I"  lle^ll 
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therefore  the  parameter  must  equal  the  spectral  radius  for 
the  approximation  to  parallel  the  error  norm  in  the  limit 
as  iteration  increases. 

Fi.qure  21  shows  a  comparison  of  the  three  error  limit 
methods  and  the  error  norm.  The  parameter  '.ms  chosen  to 
equal  the  SOU  spectral  radius.  note  that  method  three  does 
not  approximate  the  error  norm  as  well  as  methods  one  and 
two .  This  was  typical.  fee  Figures  30  and  39. 

Figure  20  shows  that  the  displacement  norm  ratio 
appears  to  oscillate  about  a  value  close  to  the  SOh  spectral 
radius.  From  Equation  (12),  the  limiting  value  of  the 
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displacement  norm  ratio  is  in  fact  equal  to  the  spectral 
radius.  Averaging  the  displacement  ratio  over  a  number  of 
previous  iterations  is  one  way  to  approximate  the  spectral 
radius . 

One  advantage  of  method  two  over  method  one  might  be 
the  fact  that  the  domain  of  the  parameter  K  for  method  two 
(Figure  5)  is  the  same  as  the  range  of  possible  convergent 
spectral  radius  values,  namely,  zero  to  unity. 

With  respect  to  additional  computer  operations,  methods 
two  and  three  require  only  the  limit  method  equation  be 
algebraically  solved.  method  one  requires  an  additional 
iteration  step  to  determine  the  nth  displacement. 

xw  =  x''’'0  +  J Cn'° 

ll&^H  *  II  dMH /O- k) 

•where  as  the  displacement  needed  for  method  two  is  already 
known 

x<">  =  f  Jfr-O 

UeMll  *  IIJ^-’ll  K/(l-K) 

and  method  three  uses  only  the  initial  displacement. 

=  x'”'0  .  J*"-0 

II e. <n> II  -  I J(0>ll  KV6-K) 

method  two  appears  to  be  the  best  choice  with  respect  to 
required  computer  operations,  accuracy,  and  the  intrinsic 


domain  of  its  parameter  being  zero  to  unity 


The  common  method  of  testing  for  convergence  is  to 
use  the  displacement  as  an  approximation  of  the  error 
between  matric  and  iterative  solutions.  The  displacement 
will  be  within  an  order  of  magnitude  of  the  error  as  long 
as  the  actual  parameter  K  associated  v/ith  method  two 
(Figure  17)  remains  between  0.1  and  0.9  (Figure  18). 
Remembering  that  K  approaches  the  spectral  radius  as  the 
number  of  iterations  become  large  suggests  that  the  dis¬ 
placement  can  be  used  when  the  spectral  radius  is  between 
0 . 1  and  0.9. 

Errors  can  be  expressed  in  absolute  or  relative  terms. 
One  is  not  intrinsically  more  a.ccurate  than  the  other.  The 
choice  is  a  matter  of  convenience  in  expression.  A  method 
which  yields  a  more  accurate  absolute  error  approximation 
v;ill  yield  a  more  accurate  relative  error  approximation. 

So  this  discussion  is  applicable  to  both  expressions. 


V  Conclusions  and  Recommendations 


Conclusions 

1.  Of  those  methods  investigated,  and  depending  on  the 
parameter  K  chosen, 

De^’l  —  K/O-K) 

represents  the  optimum  error  approximation  method  with 
respect  to  accuracy,  ease  of  computing,  and  required 
computer  resources. 

2.  The  optimum  choice  for  parameter  K  is  the  spectral 
radius  of  the  iterative  matrix.  It  can  be  calculated  using 
the  power  method  and  equations  relating  Jacobi  spectral 
radius  to  SOR  spectral  radius.  Alternately  the  spectral 
radius  can  be  approximated  by  averaging  displacement  norm 
ratios  of  previous  iterations. 

3.  The  displacement  norm  will  normally  approximate  the 
error  norm  within  an  order  of  magnitude  when  the  spectral 
radius  is  between  0.1  and  0.9. 

4.  The  spectral  radius  of  the  iterative  matrix  increases 
as  spacial  nodal  density  increases,  thus  convergence  r'ate 
decreases . 

5.  The  spectral  radius  of  the  iterative  matrix  increases 
as  the  time  interval  between  solutions  of  the  implicit 
matric  equation  is  increased.  Thus  convergence  rate 
decreases . 

Recommendations 

Since  the  Jacobi  iteration  matrix  for  transient  heat 
conduction  is  simple  and  retains  its  simple  form  as  spacial 
nodal  density  is  varied,  a  particularly  simple  similarity 
transformation  may  exist  that  diagonalizes  the  Jacobi  matrix. 
The  Jacobi  eigenvalues  and  spectral  radius  could  then  be 
obtained  from  it  rather  than  from  a  lengthy  iterative 
procedure  such  as  the  power  method. 
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Appendix  A 

Eigenvalue  Error  Limit 


Derivation  of  the  error  limit  between  the  quotient 
approximation  and  the  eigenvalue  closest  to  it  (10:462-463). 

C  .  =  orthonormal  eigenvectors  of  matrix  A. 

X  =  (c;£j)  ,  expanded  in  terms  of  eigen¬ 

vectors. 

y  =  Ax 
=  A  ?, 

I  “  I 

=  1  (C;A  £;  ) 

itl  '  ' 

Ae;  =  A  e.; 

=  I  (C;  >,•£;) 

V=  arbitrary  vector 

-  yT y  /  (v1*) 

n  O 

=  Z  (C;[X;-q  ]e.  ) 

i-i  r>  ' 

(y-%x)T(y-i'X)*  |,  (c.‘f  -<},]  e,‘) 

=  real 

C  A;  ]  =  |  A;'  ^  | 
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(y -,,x)T(y-?,*)=  I, 

Xr- A  .  closest  to  • 

|  A j.  -  ^,1  *  I  X;  -  «J.,  I  f°ra11  1 

if  Cc;‘|A^-  %,r ) 

=  ia  >-*,r  i/c,.v) 

I  ('c‘e,‘)  =  x7x 


h  >09, 


'I  A  %-<^y  m. 

IAV-^,I  -=e. 


=  e, 


(y  -  X  )  =  y  -  * ) 


=  scalar 


T  ^  T 

=  y  -  % * 


(y-^)T(r^=(y,i.x^r^x) 

=  /ryV-vV?>T' 

=  <^xTx)  -  <j/yT*  txTy)  *  Vx) 

yTX  =  symmetric,  scalar 

yTx  =  CyTx)T 


(1:162,165) 


y  Tx  «  X  Ty 
(y-<%  *y  (  y~^#x)=  (xrx)-Z<i'  (xr y)  *  (yTy) 

m,  =  x  Tx 


m,  ^  x  y 

mt=  yry 

*  -  2  ^|/>9I 

=  /'>?.  -  2  xn./sno  +  ^9t  /'jo.  ) 

~  2  ^  t  xy?x/mo) 

=™.  s; 

e,z m.  -  (  y  -  (y  -%tx)  ^  $>, 

me  >  o 

e ;  *  s,1 


Similarly 


=  orthonormal  eigenvectors  of  matrix  A. 


X  = 


£ 
;=  t 


(c:c:) 


expanded  in 


terms  of  eigenvectors. 
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I— Mftfti I ^I^iifir , nifriit if* i w 

MBHIid— mum-  i 


A1  2  (c;&;) 

I  ( C;  A'e;  ) 

»*-• 

Ac;  =  X.e; 
f  ^C;A;*C;) 

’  =  arbitrary  vector 

VT2  /  (vTx) 

Z-<j >  i  Af*C;)  -  ^ 

=  1  (c.[  A,*- 

(h  -  k)t  (i  -  "  ?,  ^c;  c’: 

v  i 

Aj  -  =  real 


Xo  =  X.  Closest  to 

>  1 


I  v  t 


I  I  \  1  _  ^ 


for  all  i 


=  l\"  r  |, 


£  (c.xc..‘)  = 


T 

X  X 


= 

=  1 V'?!1* 

i  V  ■  %i  -  ei 

-  e/m. 

x)T  =  ZT-  (<^x)T 


^  =  scalar 

_  T  T 

=  Z  -  <J£  X 

(i-  <^x)T  <* -<*,x)*<'aT-^xTXz- <^£x 

T  T  T  Z  T 

=  2  Z  -  Z  *■),»  X 

-  <f.x  (  x  Tx)  -  <^!('ztx»xtz)+  ZTZ 

T 

2  X  =  symmetric,  scalar 

2TX  =  T  x)T 

r 

-  X  2  (1:162, 

-  (xrx)  -  Z  c^^(xTz)  +  (z  T?) 
m0  =  xrx 

rr)3  *  xT? 
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*  '  2  mi  "  ry><- 

--  m,  (%  -  2  t™.  '  /n°‘  ^ 

$t  5  2  /*,  j™.  *  m./™. 

-  m.  C 

e;m.  *  (z-<^  x)TCi  -y)  =  S 


rr)0  >  o 

«,**  C 

e,  -  S, 

=  Z^'  -  2  ^  /r»,  /m,  *  /rry.  ) 

Must  now  find  a  limit  on  £3  I  -/^.)l 

s. 


from 


^  *  e,  -  U/  ~  %J  -° 

if  xf'  M*  ~  ° 

then  I  At  ~  %i  ^  ^ 


$.  * 


-  <Ll 


z 


i  o 


J(‘ j.£  <•  st)  -  J(i^)  -  *  ° 

J(%  *  V  -  J(%)  - 1\- 

£> 1  I  ><x  '  s(%)i 

J(< S.)  -  j(%j  -£iio 

s.  e2  s  i  V-  V  *° 

If  «». s  V  *  ° 

then  |  Xj1  "  ^  I  =  “}*  ~  X^  *  ° 

S.  *  <?,  -  X^1  *  O 

V  -  <*,-  st 
-  V  - 

If  ‘h  4 


then 


■*,  -  V*  <*. '  s-  so 

%>  *  k  4-0-° 

o  *  x%  -  y^t)  -  s,)-  -A^,) 

o  6  -  xt  <■  /(^.,)  -  " O  *■ 
o  —  I  xv-  /(%)!*  -'/^ k'O  * 

€,  =  I  x»  '  I 
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Appendix  B 

Ilatric  llultinlication  of  Principal  Vector's 

Multiplication  of  a  principal  vector  by  an  integer 
power  of  its  matrix.  Assume  the  rath  eigenvalue  is  multiple 
and  the  mth  eigenvector  must  be  supplemented  by  principal 
vectors  in  order  for  the  vector  space  to  be  spanned 

•  •  •  P  =  ’f’  "f*  •  *  #  P  >«  » 

I  '  nl  )  1  1*7*1  t  )  ’*n*i  )  '  m  r  jj  *•  I  7 

where 


C  ;  =  eigenvector 

t;  ,  i  >  m  =  principal  vector 

5  =  number  of  supplemental  principal  vectors 

=  ~t~ ^  =  eigenvector  from  which 

derived 

From  the  defining  relation  of  a  principal  vector  (5:30) 

At.  =1.  -•  t; 

A*t;  -  A  (At;) 

*  A  (t;.,  *■  t;) 

A  +  Am  A  t; 

-  ( A. +  A  t;_, )  *  Am  <>._,♦  t;) 

=  A,  *  2  A  A,  +  >C  t, 
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For  the 


'..'here 


A't;  =  A  (A' it) 

■A  K  -  C  t.) 

=  A  t;.,  '  2  A_  A  f,..,  -■  A.,1  A  t; 

--  (  Z-,  <-  K  h-J  *■  6,..,  *  f.J 

"  (ft,,,  ♦  A,  t.) 

■  '3A,  t.,t  '3A;t.(  -  A„J  t. 

A‘t;  =  A  ('A’O 

=  A  ^  t;.,  *  3  A„  t,..„  -  3  A,'  f;.(  <•  A^  t. ) 
■At..,  *-  3  X„  A  r, i  t3Aj  fit.,,  *A.sAt; 

f.J'JAjt,.,  *A,t..J 
^  3  A_*  if  f;.,  *  A_  t;„)  *  A„5  ( t;.(  <•  A„  r,) 

=  t  ■  +  4-  A  t  •  +-  6  x  t . 

1-4-  'n  <  -  z. 

m3  t.  *  X„T  t. 

i  —  /  / 


ICth  pov.'er  of  matrix  A 


A*  t;  -  Z  L  (j)  \„*-J 


A  =  matrix 

K  =  the  power  to  v/hich  the  matrix  is  raised 
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=  principal  vector, 


l  ^  no 


t. 

i 


j  =  0/  l-m  if  K  -  l -nn 
j  =  O,  K  if  K  -  l-  no 

=  e. ^  =  eigenvector  from  which  t  - 

is  derived 


(j)  ~  (  K)  j  /(  J)  \/i^~  j)  f  =  binomial  coefficient 
(3:624) 


Appendi 


Solve 
(4:32-35)  . 


Append i;:  D 

Validation  of  Analytic  Solution 
one  dimensional  heat  conduction  equation  first. 

i =  (i/k)  3t /it 

T  =  AW6W 

yr/ix‘  *  a"(x)  s(t) 

)T /it  =  A(x)  6 VO 

A' A)  6(0  =  (</*)  A (x)B'(t) 

\' (*)/*(*■)*( i / k)  =  -x8- 

A’W'-k'aW 

6  VO  r  -  *  *  k  6  A) 

bVO  =  -  l1  B (t) 

8  A)  -  c. 

6 VO  =  -  Se'St 

c  -  frf  .  «■  -  &  f 

-  a  e  =  -  A  e 

S  =  A* 

6(tJ=e-J'' 
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*  =  n 


T0  f  *,  fO  =  ^e^-c'^)cae-A  1 


-  C4  sin  (nx)  e. 


-A*t 


Apply  initial  conditions 

S/o  (x  )  =  Tn  (x,  o) 
-  C+  5/0  (nx) 


1 


>0  ^  t 


"T (x,  t)  ~  sun  (x)  e. 


T( x/  tJ  =  t>in  (x)  e- 


-/k  f 


<*>  -  n 


T tJ  -  5<o  (*)  e. 


n  =  1 


TO,  t)  -  &>n  (k)  q. 


-Kt 


Solve  two  dimensional  problem  using  one  dimensional 
solutions  (2:33-34). 


T  (*,y,  f)  *  T  T  Sy,  t) 

-Z  Kt 

=  S  m  (k)  S"?  Cy)  e 

Boundary  conditions  are 

T  (0,y,t)  -T  (n,  y,r) 

-  T (x,0,t)*  T(*,ir,t)  =  o 

Initial  conditions  are 

T  ^  X,  y,  <?)  ~  sin  (*)  sinCy) 

Doiiiain 

o  -  x  -  7T 
O  6  y  6  TT 

0±  t 


Validation  of  solution 


5m  CK)  s>n(y)  C- 
}T /^X  *  5 >>n(y)  O 

-5 in  (x)  sin  Cy)  & 


/}ylz  -  sin  (*)  sin  Cy)  e. 


—  2.  K  T 


znt 


-  2« 


-ZKt 


-  2  Kf 


)T/)t  -  -  ZK  sin(x)  S*n(y)  & 

yi/^y7'  =  sin(x)sinCy) 


-  2  K 

& 
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Appendix  F 

Computer  Code-  Spectral  Radius,  Double  Iteration 
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Appendix  N 


Iterative  Error  Approximations 
3  by  3  Nodal  Density 

Number  of  interior  nodes  =  4 
Time  interval  between  matric  solutions  =  1 
Diffusion  coefficient  =  1 
Computer  precision  =  27  decimal  dibits 
Number  of  spectral  radius  iterations  =  30 
Jacobi  spectral  radius  =  .33242 
Gauss-Seidel  spectral  radius  =  .15399 
SOR  optimum  spectral  radius  =  .041702 
Number  of  iterations  =  30 


?uro  23.  i’araneter  Required  for  Error  Li, -nit  ilethod  1  to  Equal  the  Error  iiorm 
3  x  3  i Social  Density 


LEGEND  (TOP  TO  BOTTOI 


24.  Order  of  iiar;nitude  Comparison  between  Error  Norm  and  Displacement  Norm 
Error  Limit  iiethod  1,3x3  Nodal  Density 


Figure  25.  Error  norm  to  Displacement  norm  Ratio,  Error  Limit  Method 
3x3  Ilodal  Density 


( I+N3/N0) /I -M 


ure  26.  Parameter  Required  for  Error  Limit  method  2  to  Equal  the  Error  Horn 
3x3  Ilodal  Density 


Order  of  Magnitude  Comparison  between  Error  Norm  ancl  Displacement  Norm 
Error  Limit  iiethod  2,  3  x  3  Nodal  Density 


ure  28.  Comparison  of  Error  Norm  to  Error  Limit  Method  3,  3  ;;  3  Nodal  Density 


Figure  20.  Displacement  Norm  Ratio,  3x3  Nodal  Density 


-rure  30.  Three  Error  Limit  Methods  Compared  to  Error  Norm,  E  =  SOU  Optimum 
Spectral  Radius,  3  ;;  3  Modal  Density 


Appendix  I 

Iterative  Error  Approximations 
20  by  20  Modal  Density 


Humber  of  intex'ior  nodes  =  361 

Time  interval  betv/een  matric  solutions  =  10,000,000 

Diffusion  coefficient  =  1 

Computer  precision  =  13  decimal  dibits 

Humber  of  spectral  radius  iterations  =  30 

Jacobi  spectral  radius  =  .93715 

Gauss-Seidel  spectral  radius  =.97447 

S0R  optimum  spectral  radius  =  .72446 

Number  of  iterations  =  200 
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ement  Norm  Ratio,  Error  Limit  method 


Order  of  Magnitude  Comparison  between  Error  Norm  and  Displacement  Norm 
Error  Limit  Method  1,  20  x  20  Nodal  Density 


ir.ure  35.  Parameter  Required  for  Error  Limit  Method  2  to  Equal  the  Error  Norm 
20  x  20  Nodal  Density 
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ilorm  to  Error  Limit  Method 


•-■•*** 


Figure  38.  Displacement  tloru 


Figure  39.  Three  Error  Limit  llet’nods  Compared  to  Error  Norm,  K  =  SOU  Optimum 
Spectral  Radius,  20  ;;  20  Nodal  Density 
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